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Abstract 

The spatial variability of stress fields resulting from polycrystalline aggregate calculations in¬ 
volving random grain geometry and crystal orientations is investigated. A periodogram-based 
method is proposed to identify the properties of homogeneous Gaussian random fields (power 
spectral density and related covariance structure). Based on a set of Hnite element poly crystalline 
aggregate calculations the properties of the maximal principal stress field are identified. Two cases 
are considered, using either a fixed or random grain geometry. The stability of the method w.r.t 
the number of samples and the load level (up to 3.5 % macroscopic deformation) is investigated. 

Keywords: Polycrystalline aggregates - Crystal plasticity - Random fields - Spatial variability 
- Correlation structure 


1 Introduction 


In pressurized water reactors of nuclear plants, the pressure vessel constitutes one element of the second 
safety barrier between the radioactive fuel rods and the external environment. It is made of 16MND5 
(A508) steel which is forged and welded. In case of operating accidents such as LOCA {loss of coolant 
accident)^ the pressure vessel is subjected to a pressurized thermal shock due to fast injection of cold 
water into the primary circuit. If some defects {e.g cracks) were present in the vessel wall this may 
lead to crack initiation and propagation and the brittle fracture of the vessel. The detailed study of 
the embrittlement of 16MND5 steel under irradiation is thus a great concern for electrical companies 
such as EDF. 

The brittle fracture behavior of the 16MND5 steel has been thoroughly studied in the last decade 


using the local approach of fracture theory (Tanguy 2001) and the so-called Beremin model (Beremin 


19831, which assumes that cleavage is controlled by the propagation of the weakest link between a 
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population of pre-existing micro-defects in the material. This approach has been recently coupled with 
polycrystalline aggregates simulations (Mathieu et al, 20061, ( Mathieu et al[ 20101. 

The main idea is to model a material representative volume element (RVE) as a polycrystalline 
synthetic aggregate and compute the stress field under given load conditions. As a post-processing a 
statistical distribution of defects (carbides) is sampled over the volume. In each Gauss point of the 
finite element mesh the cleavage criterion is attained somewhere along the load path if a) the equivalent 
plastic strain has attained some threshold (cleavage initiation) and b) a GrifRth-like criterion applied 
to the size of the carbide in this Gauss point is reached (cleavage propagation). Within the weakest 
link theory the failure of a single critical carbide induces the failure of the RVE. 

From a single RVE simulation (i.e. a single stress held) various distributions of carbides are drawn, 
each realization leading to a maximal principal stress associated to failure. Then the distribution of 
these quantities is htted using a Weibull law ( Mathieu et al} 20101. In such an approach, the current 
practice of computational micromechanics assumes that the RVE is large enough to represent the 
behavior of the material so that a single polycrystalline analysis is carried out (the large CPU required 
by polycrystalline simulations also favours the use of a single simulation). However it is believed that 
numerous parameters such as grain geometry and orientation may inhuence the stress held and thus 
the hnal result. 

The connection between micromechanics and stochastic methods has been given much attention 
in the past few years, as shown in Graham-Brady et al (2006); Stefanou (2009); Xu and Chen (2009). 
Many papers are devoted to developing probabilistic models for reproducing a random microstructure, 


e.g. 

Graham and Baxter 

(2001); 

Liu et al 

(2007) 

Ghung et al ( 

2005 

); Chakraborty and Rahman 

(2008 

). The specihc representation of polycrystalline microstructures has been addressed in 

Arwade 

and Grigoriu 

(2004) 

Grigoriu 

(2010|; Li et al 

(2010) 

Kouchmeshky and Zabaras 

(2010 

) among others. 


The propagation of the uncertainty on the microstructure through a micromechanical model in order 
to study the variability of the resulting strain and stresses has not been addressed much though (see 

e.g. 


Kouchmeshky and Zabaras (2009)). 


In this paper it is proposed to identify the properties of a stress random held resulting from the 
progressive loading of a polycrystalline aggregate. More precisely, assuming that the stress random 
held is Gaussian, a procedure called periodogram method is devised, which allows one to identify the 
correlation structure of the resulting stress held. 

The paper is organized as follows: in Section 2 basics of Gaussian random helds are recalled and the 
periodogram method is presented (Dang et al 2011b). The polycrystalline aggregate computational 
model is detailed in Section 3. The methodology for identifying the correlation structure of the resulting 
stress held is presented in Section 4. Two application cases are then investigated, namely an aggregate 
with hxed grain boundaries and random crystallographic orientations (Section 5) and an aggregate 
with both random geometry and orientations (Section 6). The variance of the resulting stress held 
as well as the spatial covariance function and its correlation lengths is investigated in details. The 
properties of the identihed random helds will be used in a forthcoming study in the context of the 
local approach to fracture, as explained above. 


2 Inference of the properties of a Gaussian random field 


In this section an identihcation method called periodogram is presented, which uses an estimator 
of the Power Spectral Density (PSD) in order to identify the correlation structure of a Gaussian 


homogeneous random field. Based on original developments by 

Stoica and Moses 

19971 

and 

Li 

(2005 

for unidimensional helds, it has been extended to two-dimensional cases by 

Dang et al| ( 

2011b) 

. As it 


relies upon the use of the Fast Fourier Transform (EFT) its computational efficiency is remarkable. 
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2.1 Definitions 


A Gaussian random field Z{x,uj) is completely defined by its mean value its standard deviation 

cr(a;) and its auto-covariance function C{x,x'). It is said homogeneous if the mean value fJ.{x) and the 
standard deviation a{x) are constant in the domain of definition of x and the auto-covariance function 
C{x,x') only depends on the shift h = x — x'. Let us introduce the n—th statistical moment and 
the spatial average 


OO 

=E[Z”'{xo,uj)]= j z'^{xo,uj)fziz,Xo)dz 

— OO 

= lim ^ [ z"{x,ujQ)dx 
V— >-oo V J 


( 1 ) 

( 2 ) 


The field is said ergodic if its ensemble statistics is equal to the spatial average, i.e. 


(Cramer and Leadbetter 1967). Several popular covariance models for two-dimensional homogeneous 


random fields are presented in Table In this table, a is the constant standard deviation of the field, 
hi, h 2 are the components of the shift h in the two directions, li, I 2 are the correlation lengths in the 
two directions. Gaussian and exponential models are plotted in Figure for the sake of illustration. 
Note that we call correlation length the parameters that appear in the definition of the covariance 
functions. This is not to be confused with the scale of fluctuation ( [Vanmarcke 1983), which combine 
both the shape of the covariance function and the lengths fi, l 2 - In one dimension, denoting by p(x; 1) 
the autocorrelation function, the scale of fluctuation may be defined by: 


2L = 


p{x; 1) dx 


which reduces to Ic = I for the exponential correlation function and Ic = y/^l2 ~ 0.886 Z for the 
Gaussian case. Similar expression are available in two and three dimensions, see e.g. |Xu and Chen| 
(l2009|. 


Table 1: Govariance functions and associated power spectral densities for homogeneous twodimensional 
random fields 


Model Govariance function Power spectral density 

Exponential cr^exp 

Gaussian cr^exp 

Wave (T^sinc 

Triangle (T^tri(f 

_ Wi ^ /2 

(hi _|_ ^2 A 

V 72 “T ,2 ) 

^)tri(M) 

-) 

) 

2li 2 I 2 

^ T+4F2]f7f l+47r2i^/| 

TrZiexp 7rZ2exp (tt^ZI/I) 

cr^ 7rZirecti(7rZi/i) 7rZ2recti(7rZ2/2) 

Zisinc^(7r/iZi) hsinc^{ tt f2l2) 


sinc(x) = s\a.x/x 

tri(x) = l—\x\ if |a;| < 1 and 0 otherwise 
lectrif) = 1 if I/I < i and 0 otherwise 


The power spectral density (PSD) of the random field is the Fourier transform of its covariance func¬ 
tion as a result of the Wiener-Khintchine relationship (Preumont 19901. The following relationships 
hold: 


S{fi,f2)= f f 

— OO —OO 

OO OO 

Cihi,h2)= J / 

— OO —OO 

The PSD of the Gaussian and exponential covariance models are presented in Table 


(3) 

(4) 
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Figure 1: Gaussian covariance model (left) and Exponential covariance model (right) : a = 2, 
li = I 2 = 5 


2.2 Identification of a PSD 

2.2.1 Empirical periodogram 

One considers an ergodic homogeneous random field Z{xi,X 2 ), {xi,X 2 ) G x I ?2 C for which a 
single realization z{xi^X 2 ) is available. If the random field was defined over an infinite domain, the 
classical estimation of the covariance function would be: 

^ 7V-1 M-1 

Z{xin + hi,X2m + h2)Z{xin,X2m) (5) 

n— — N m——M 


By definition, the Fourier transform of the covariance estimation is an estimation of the PSD. 

5 (/ i ,/ 2 ) = ^Zih,h)Z*if,,h) = ^\ZihJ,)\^ ( 6 ) 

where |.| denotes the modulus operator. 

In practice, the problem is to estimate the periodogram from a limited amount of data gathered on 
N X M grid {z{xii, X 2 j), i = 1, ... ,N;j = 1, ... , M}. Due to symmetry, the covariance estimation 
in Eq.([^ is recast as follows: 

N-kM-l 

= WM S S Z{xin+hik,X2Tn + h2l)Z{xin,X2m) (7) 

n—1 m—1 


By taking the expectation of the above equation one gets: 


E 


CQiik, h2i) 


N-kM-l 

~/V 


E \^Z(Xin “f l^lk} X2m 


N-kM-l 

~/V 


C[h^,h2) 


Z i^Xin: 3:2m)] 


( 8 ) 


The latter equation exhibits some bias term between the expectation of the estimator and the covari¬ 
ance function C{hi, h 2 )- Using the symmetry of the covariance function, one can write: 


E 


C{hik, h2i) 


WB{k,l)C{hi,h2) 


(9) 
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where wsik, 1) is the triangle window, also known as the Bartlett window (Figure]^: 

if \k\ < N; |1| < M 


WB{k,l) = 1 Q ^ 


M 


otherwise 


Consequently the expectation of the periodogram estimation becomes: 


E 


S{flk,f2l) 


= T 


{e [C{hik, h2i)] } = HwB{k, l)C{hi,h2)} 
= WB{flj2)*Sihj2) 


( 10 ) 


( 11 ) 


where if and VFb(/i, / 2 ) respectively denote the 2D Fourier transform operator and the Fourier trans¬ 
form of the Bartlett window and * denotes the convolution product. This window tends to a Dirac 
pulse when N, M tend to infinity and wb tends to a unit constant. Thus the periodogram estima¬ 
tion is asymptotically unbiased. However it is not consistent since its variance does not tend to zero 


(Preumont 1990). Furthermore using this window leads to a convolution product which introduces 


additional computational burden. Hence in practice, the modified periodogram presented in the next 
section is used to estimate the PSD of the random field. 


2.2.2 Modified periodogram 

The modified periodogram is built up in order to avoid the convolution product with the transformed 
window VFb(/i,/ 2 ) in Eq. In this approach, the data is multiplied directly with the window 

w{x,y) before the Fourier transform is carried out. It aims at filtering the data to limit the influence 
of long distance terms and to focus on the information given by the short distance terms. This leads 
to the following estimate: 

S{fl,f2) = {z{xi,X2).w{xi,X2)} P (12) 

where U is the energy of the window calculated by: 


U = 


1 


N M 


Di D 2 


2^2^W {xu,X2j) 


(13) 


i=l j=l 

and Di^D 2 denote the size of the two-dimensional domain Vi x I? 2 - Various window functions are 
proposed in [Preumont (1990), see Table In this paper we will use mainly the Blackman window 
(Figure [ 2 ]). 

2.2.3 Average modified periodogram 


As shown in Section 2.2.1 the estimation of the periodogram is asymptotically unbiased, however not 
consistent since its variance does not tend to zero when N, M tend to infinity. The averaging of the 
modified periodogram will solve this problem. Assume that L realizations of th e ra ndom field are 
available. For each realization z''{xi,X 2 ), one calculates the periodogram as in Eq.(12): 

S\fi,f2) = 1-^ {z\xi,X2).w{xi,X2)} P (14) 

with 1 < I < L. Then one calculates the average periodogram by: 

1 ^ 

5(/i,/2) = yE^'(/i’/2) 


(15) 




Therefore the variance of the average periodogram is: 

Varp(/i,/ 2 )] = ^Var[5(/i,/2) 


(16) 


It is then obvious that this variance tends to zero when L tends to infinity, making the “average 
modified periodogram” approach more robust. 
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Model 


Table 2: Window functions used in the modified periodogram approach 
Window equation 


Bartlett 


Hann 


Hamming 


Blackman 


/ if |fc| < iV; |/| < M 

\ 0 otherwise 

r [0.5 + 0.5cos(^)] [0.5 + 0.5cos(||)] if |fc| < N; \l\ < M 
\ 0 otherwise 

r [0.54 + 0.46cos(^)] [0.54 + 0.46cos(||)] if |fc| < N- |Z| < M 
\ 0 otherwise 

f [0.42 + 0.5cos(^) + 0.08cos(^)] [0.42 + 0.5cos(^) + 0.08cos(^)] 

I 0 


if |fc| < N; \l\ < M 
otherwise 



Figure 2: (a) Bartlett window ; (b) Blackman window 


2.2.4 Final algorithm for PSD estimation 

As a summary, the algorithm to estimate the PSD of a random field from L realizations may be 
decomposed into the four following steps: 

1. multiplication of each realization by a selected window, e.g. the Blackman window (see Table[^; 

2. computation of 2D Fourier transform of the product of the current realization by the filtering 
window; 

3. computation of the modulus of the result to obtain the PSD estimation of each realization; 

4. averaging of the L PSD estimations. 

Once the empirical periodogram has been computed, a theoretical periodogram is selected {e.g. 
Gaussian, exponential, etc., see Tableand fitted using a least-square procedure (Marquardt, 1963). 
In case of multiple potential forms for the theoretical periodogram the best fitting is selected according 
to the smallest residual. 


(Marquardt, 1963 
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3 Modeling polycrystalline aggregates 

In this section the computational mechanical model used in this study is presented. It simulates a 
tensile test on a bidimensional polycrystalline aggregate under plane strain conditions. The various 
ingredients are discussed, namely: 

• the microstructure of the material and its synthetic representation; 

• the material constitutive law; 

• the boundary conditions applied onto the aggregate; 

• the mesh used in the finite element simulation. 


3.1 Material characterization 


The material is a 16MND5 ferritic steel with a granular microstructure. The ferrite has a body centered 
cubic (BCC) structure. Three families of slip system should be taken into account, namely {1I0}(1I1), 
{112}(111), {I23}(1I1). However, following Franciosi (1985) it is assumed that the glides on the plane 
123 are a succession of micro-glides on the planes 110, 112. This leads to consider only the two first 
families, which yields 24 slip systems by symmetry. 


3.2 Crystal plasticity 


The model for crystal plasticity chosen in this work has been originally formulated in Meric and 


Cailletaud (19911 within the small strain framework. The total strain rate is classically decomposed 


as the sum of the elastic strain rate and plastic strain rate ef 


.. = o'?. -I- 


(17) 


The elastic part follows the Hooke’s law and the plastic part is calculated from the shear strain rates 
of the 24 active slip systems. 

24 


e 


p 




(18) 


where 7 ® is the shear strain rate of the slip system g and i?® is the Schmid factor which presents the 
geometrical projection tensor. The latter is calculated from the normal vector to the gliding plane n 
and the direction of gliding m. 


= -{miTij + mjTii) 


(19) 


The Resolved Shear Stress (RSS) r® of the slip system g is the projection of the stress tensor via the 
Schmid factor. 

( 20 ) 

The shear strain rates 7 ® of each slip system g are the internal variable that describes plasticity. The 
evolution of these variables depends on the difference between the RSS and the actual critical RSS 
r® in an elastoviscoplastic setting: 


/t- 9 - r9\” 

where K and n are material constants, and sign(a) = a/|a| if a 7 ^: 0 and 0 otherwise. Note that this 
formula corresponds to an elastoviscoplastic constitutive law but the viscous effect will be negligible 
if sufficiently large values of K and n are selected. Its power form allows one to automatically detect 
the active slip systems. All the systems are considered active but the slip rate is significant only if the 
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RSS r® is much higher than the critical RSS r®. This procedure allows one to numerically smooth the 
elastoplastic constitutive law. 

The critical RSS r® evolves according to the following isotropic hardening law: 

24 

r ! = + Q ® E ^"'(1 - 6 -''“’"—) ( 22 ) 

S = 1 


where 7 *^^ = I The exponential term presents the hardening saturation in the material when 

to 

the accumulated slip is high. r®Q is the initial critical RSS on the considered system g. Q® and 6® 
are parameters which depend on the material, /i®® is the hardening matrix of size 24 x 24 whose 
component ft,®® presents the hardening effect of the system g on the system s. In the present work, one 
considers only two families of slip systems named 110(111), 112(111). Thus the hardening matrix ft®® 
is completely defined by four coefficients fti,ft 2 ,ft 3 ,ft 4 only. The values of these coefficients and this 
matrix are presented in Mathieu (2006). All the parameters describing crystal plasticity for 16MND5 
steel are gathered in Table 


Table 3: Parameters of the crystal plasticity constitutive law for the 16MND5 steel (Mathieu 2006) 


Isotropic elasticity Viscoplasticity Isotropic hardening 


E 

{MPa) 

V 

K 

{MPa.s^/'^) 

n Tea 

{MPa) 

Q b hi h2 fta 

{MPa) 

210,000 

0.3 

15 

12 175 

20 30 1 1 1 1 


3.3 Microstructure and boundary conditions 


The construction of the aggregate is based on the Voronoi polyhedra model (Gilbert 1962), generated 


in this work with the Quickhull algorithm (Barber et al 1996 1 . The geometry of the resulting synthetic 


aggregate, which is a simplified representation of the real microstructure of the 16MND5 steel, is shown 
in Figure]^ It corresponds to a square of size 1,000 (this is a relative length which shall be scaled with 
a real length depending on the grain size). Grain boundaries are considered as perfect interfaces. Note 
that more detailed models of grains have been proposed recently using so-called Laguerre tessellations 


(Lautensack and Zuyev 

2008 

1 in order to better fit the observed distributions of grain size, see e.g. 

Zhang et al. ( 

20111 ; 

Leonardia et al. 

(2012 

1 . 


The same crystallographic orientation, defined by the three Euler angles ipi, 0, :/? 2 , is randomly 
assigned to all integration points inside each individual grain using a uniform distribution. In Figure 
[^a, the color of each grain corresponds to a given crystallographic orientation. The mesh is generated 
by the BLSURF algorithm (Laug and Borouchaki 1999) of the Salome software (http://www.salome- 
platform.org). The mesh of the generated specimen is presented in Figurej^b. The finite elements are 
quadratic 6 -node triangles with 3 integration points. 

The boundary conditions applied onto the aggregate are sketched in Figure]^ The lower surface 
is blocked along the Y direction. The displacements DX = DY = 0 are blocked at the origin of the 
coordinate system (lower left corner). On the upper surface, an homogeneous displacement is applied 
by steps in the Y direction up to a macroscopic strain equal to 3.5%. The computation is carried out 
using the open source finite efement software Code_Aster (http://www.code-aster.org). 

The computational cost for such a non linear analysis is high. The number of degrees of freedom of 
the finite element model is 33, 885. A parallel computing method based on sub-domain decomposition 
is used. One simulation of a full tensile test up to 3.5% strain requires about 2 hours computation 
time when distributed over 4 processors. 





































Figure 3: (a) Two-dimensional polycrystalline aggregate modelling a volume of 16MND5 steel (100 

grains) (b) Mesh of the specimen (11,295 nodes) 


Y 



DY - 0 


Figure 4: Boundary conditions used for simulating the tensile test 


3.4 Results 


In this section, we present the result of the simulation of a tensile test on the 2D aggregate at different 
scales. We define the mean stress and strain tensor calculated in a volume V by: 


^ = ^ J 

(23) 

V 


E=i|cciF 

(24) 


V 
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Figure 5: Macroscopic strain-stress relationship in the X and Y directions 


Figure shows the macroscopic strain/stress curve. It is observed that Sxx = 0 as expected 
whereas the uniaxial behaviour shows a global elastoplastic behaviour. 

At the mesoscopic scale one can observe the mean strain-stress relationship in each grain as shown 
in Figure Because of the different crystallographic orientations in each grain, the mean elastoplastic 
behaviour is different from grain to grain. Furthermore, whereas the mean stress Yxx calculated in 
all the specimen is zero, the mean values calculated in each single grain are scattered around zero. 
This observation shows the first scale of heterogeneity of the material. 



Figure 6: Mesoscopic behavior in each grain in the X (left) and Y (right) directions 
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The microscopic behaviour of a single grain (Grain #24, see tag in Figure]^ is finally studied. The 
mean behavior and the strain-stress relationship at each node of this grain are plotted in Figure]^ for 
four levels of macroscopic strain, namely Eyy = 0.15%, 0.65%, 1.5%, 3.5%. 



Figure 7: Microscopic strain-stress relationship for various nodes within Grain #24 and mean tensile 
curve 


In this figure the blue point represents the stress field within the grain for a macroscopic strain 
level Eyy = 0.15%. This single point shows that the stress field is homogeneous within the grain 
in the elastic domain. The red points represent the strees values uyy in each node of the grain at 
Eyy = 0.65% macroscopic strain. One observes that the mean strain calculated for this single grain is 
0.85% and the maximal strain value £yy in a specific node may attain about 2%. Similar effects are 
observed at other levels of macroscopic strain, which show the heterogeneity of the strain and stress 
fields at the microscopic scale. It is observed that the scattering around the mean curve increases with 
the macroscopic strain. Indeed, for the final loading step corresponding to Eyy = 3.5% the mean 
strain in the grain is about 4.54%, while the local strain varies form 2.4 to 9%. 


4 Identification of the maximal principal stress field 

In this section the method developed in Section 2 is applied to the identification of the properties of the 
random stress field in polycristalline aggregate calculations. More specifically the maximal principal 
stress field aj that is computed from repeated polycrystalline simulations is considered. Throughout 
the paper this stress field is considered Gaussian. This is a strong assumption which shall be considered 
as a first approximation. Indeed the maximal principal stress is positive in nature under the uniaxial 
loading that is considered and a Gaussian assumption cannot totally fit this feature. Yet it is believed 
that the results obtained in terms of the description of the spatial variability (covariance functions), 
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which is the main outcome of the paper, will not be strongly influenced by this assumption. Note that 
methods for identifying the properties of non Gaussian random fields have been recently developed, 
see e.g. 


Perrin et al. (20141. 


4.1 Finite element calculations and projection 

The maximal principal stress field is assumed to be Gaussian and homogeneous (the latter assumption 
will be empirically checked as shown in the sequel). The periodogram method is applied using K = 
realizations of stress fields, i.e. 35 full elastoplastic analysis of aggregates up to a macroscopic strain 
of 3.5 %. The identification is carried out successively at various levels of the macroscopic strain. Two 
cases are considered: 


• Gase #1: the grains geometry is the same for all the finite element calculations. Only the 
crystallographic orientations are varying from one calculation to the other. 

• Gase ^2: both the grains geometry and the crystallographic orientations vary. 


The input data of the identification problem is the maximal principal stress field cr/ obtained from 
the finite element calculations. As the periodogram method is based on a regular sampling of the 
random field under consideration, the brute result {i.e. the maximal principal stress at the nodes of 
the mesh) has to be projected onto a regular grid. This operation is carried out using internal routines 
of Code_Aster. Note that a slice of width 100 {i.e. 10% of total size) is discarded along the edges 
of the aggregate in order to avoid the effect of boundary conditions on the computed stress field, as 
suggested in Mathieu (20061. A typical maximal principal stress held is shown in Figure]^ 




Figure 8: A realization of the maximal principal stress held cr/ 


4.2 Check of the homogeneity of the field 

As it was described in Section 2 the periodogram method assumes that the random held under con¬ 
sideration is homogeneous. From the available realizations SIGi{x,y),i = 1, ... ,35 one hrst checks 
empirically this assumption using the following approach: 

• The ensemble mean and variance of the held is computed point-by-point throughout the grid for 
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an increasing number of realizations K = 2, ^ 35: 


K 


tJ-K{x,y) = —^SIGi{x,y) 




1 ^ 

<XKix,y) = {SIGi{x, y) - y{x, y)f 


(25) 


(26) 


i=l 


If the field is homogeneous these quantities should tend to constant values that are independent 
from the position (x, y) when K tends to infinity. 

• In order to measure the magnitude of the spatial fluctuation of the latter, the spatial average 
and spatial variance of a realization of a field Z{x,y) sampled onto a N x M grid is defined by: 

N M 

1 

MZ = 


ct | = 


N M 


NM 


^=l j=i 

whereas the associated “spatial” coefficient of variation is defined by: 

CVz = ^ 

IJ-z 


(27) 

(28) 

(29) 


The spatial coefficient of variation of the ensemble mean and variance (Eqs.(25l-(26)) are com¬ 


puted and plotted as a function of K. If the underlying random field is homogeneous it is 
expected that the curves of and CV ^2 converge to zero. 


4.3 Choice of theoretical periodograms and fitting 

From a visual inspection of the obtained empirical periodograms it appears that a Gaussian or an 
exponential model of periodogram such as those presented in Table may be consistent with the data. 
However it appeared in the various analyses that the peak of the periodogram is not always in zero. 
An initial frequency is thus introduced which shifts the theoretical periodogram. Finally, due to lack 
of fitting of the single-type periodogram {e.g. Gaussian and exponential), a combination thereof is also 
fitted. The most general model finally reads: 


S{f^,fy)=ajTrGiexp Tr^lhifx - fiof /yiexp 

‘^1x2 


lllify-C) 


+ a|- 


21 


y2 


(30) 


1 + ^T^%2{fx - fxof 1 + ^^^ll2{fy - fyo’y 


(2)^ 


where lxi,^yiTlx 2 ,ly 2 are correlation lengths in each direction X and Y (aniso-tropic field) for each 
component (I) (Gaussian part) and (2) (exponential part). Similarly /ic)\/yo^/ioEyo^ are initial 
shift frequencies. 


Note that Fq.(30) corresponds only to positive values of fx,fy The periodogram is then extended 


by symmetry for negative frequencies. In terms of associated covariance models, the linear combination 
of periodograms leads to a linear combination of covariance models. The initial frequency shift in the 
periodogram leads to oscillatory cosine terms in the covariance by inverse Fourier transform: 


C{hx,hy) =(Tiexp 


V 72 /2 / 




cos{2TTf^l\x)cos{2TT6lhy) 


a^exp 


/ I \hy\ 

^ 1x2 I 


y2 J 


(31) 


cos{2TTfilhx)cos{2nfj^l'>hy) 
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In order to com pare the various fittings the least-square residual between the empirical periodogram 
S{fx,fy) (Eq.(15l) and the fitted periodogram {fx, fy) is finally computed. The following non 


dimensional error estimate is used: 


\ 


N M 


H H fyj) - V max S{fx, fy) 

i=l j=l 


(32) 


5 Results — Case ^^1: fixed grain geometry 

5.1 Check of the homogeneity 

First the homogeneity of the maximal principal stress held is checked using the methodology proposed 
in Section 4.2 
and it is seen 


Figured shows the evolution of and CV„:i . These quantities regularly decrease 




that they would tend to zero if a larger number of realizations was available. This leads 
to accepting the assumption that the random held is homogeneous since the huctuations around the 
constant spatial average tend to zero when K increases. 


Evolution of CV of the empirical mean Evolution of CV of the empirical SD field 




Figure 9: Case #1: Evolution of and CV^ with respect to the number of realizations 


5.2 Identification of periodograms at 3.5% macroscopic strain 


The average empirical periodogram obtained from L = 35 realizations of the maximal principal stress 
held CT/ at 3.5% of macroscopic strain is plotted in Figure 10-a. 


Table presents the results of the htting of the average empirical periodogram calculated from 
35 realizations of the held using three models, namely Gaussian, exponential and a mixed “Gaussian 
-I- exponential” as in Eq.([M|. 

From the results in Table it appears that the mixed model provides a signihcantly smaller least- 
square error than that obtained from the Gaussian and exponential models respectively. The corre¬ 
sponding htted periodogram is plotted in Figure [T0|b. 

In order to better appreciate the quality of the htting, two-dimensional cuts of the empirical (resp. 
htted) periodogram are given in Figures 


for two values of fy = 0; 0.0013. Figure 


11 


12 


13 Figure 11 corresponds to a cut along the X direction 


corresponds to a cut along the Y direction for two values 


of /x = 0; 0.0013. Finally Figure 13 corresponds to a cut along the diagonal fx = f 


V 


From the above hgures it appears that the htting of the empirical periodogram by a mixed model 
is remarkably accurate. It is interesting to interpret the htted parameters reported in Table First it 
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Figure 10: Case #1: (a) Average empirical periodogram of the stress field at 3.5% macroscopic strain 
- (b) best fitted periodogram 


Table 4: Fitted parameters and error estimates for the three htted models: Gaussian, exponential and 
mixed “Gaussian + exponential” 


Model 

e 





Gaussian 




Exponential 


(Eq.( 

32 

1) 

CTl 

^xl 

lyl 

f(l) 

JxO 

f(i) 

JyO 

0-2 

^x2 

ly2 

A‘2) 

JxO 

TGI 

JyO 

Gaussian 

0.0043 

69.4 

104.6 

102.9 

0.00287 

0 

_ 

_ 

_ 

_ 

_ 

Exponential 

0.0039 

_ 

_ 

_ 

_ 

_ 

84.2 

73.8 

87.5 

0.00275 

0 

Mixed 

0.0017 

54.7 

138.4 

159.1 

0.00244 

0 

57.6 

57.5 

63.5 

0.00562 

0.0028 




Figure 11: Case #1: Cut of the periodograms in the X direction 


is observed that the amplitude of each component of the mixed periodogram is similar since cti « CT 2 - 
The variance of the field is equal to af + « 6,309 MPa^. The associated standard deviation is 
79.4 MPa. As the mean principal stress is 720 MPa at 3.5% macroscopic strain, the coefficient of 
variation of the field is about 11%. 

In order to interpret the correlation length parameters let us define the mean size of a grain Sg such 
a two-dimensional aggregate. As the volume of edge length equal to 1,000 corresponds to 100 grains. 
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Figure 12: Case #1: Cut of the periodograms in the V direction 



• empirical periodogram 
— fitted periodogram 



1.2 



Figure 13: Case #1: Cut of the periodograms along the diagonal = fy 


the equivalent diameter of a single grain reads: 


D 


g ~ 



4 1,000 X 1,000 

TT 


100 


112.8 


(33) 


Thus the correlation lengths obtained from the fitting vary from 0.55 to 1.3Dg. This shows that the 
characteristic dimension of the underlying microstructure {i.e. Dg) is of the same order of magnitude 
as these parameters. In other words the scale of local fluctuation of the stress field is related to the 
grain size, as heuristically expected. Moreover, it appears that the lengths in the X and Y directions 
are almost identical. The stress field does not show any significant anisotropy in this case. 


5.3 Influence of the number of realizations 


In this section the stability of the fitted parameters as a function of the number of available realizations 
K used in the average periodogram method is considered. In practice the procedure applied in the 
previous paragraph is run using if = 8,9, ... , 35 realizations of the stress field. The evolution of the 
standard deviations (cri,cr 2 ) is shown in Figure 


14 


The evolution of the correlation lengths l(x,y){i, 2 ) 
is shown in Figure The evolution of the initial frequencies f^x 'y)o shown in Figure 
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Evolution of sigma w.r.t the data number 



Figure 14: Case #1: Evolution of the fitted standard deviations with respect to the number of real¬ 
izations K = 8, ... ,35 


Evolution of lx w.r.t the data number 



Evolution of ly w.r.t the data number 



Figure 15: Case #1: Evolution of the fitted correlation lengths in the X, Y directions with respect to 
the number of realizations K = 8, ... ,85 


From these figures it appears that the fitted parameters tend to a converged value when at least 
20 realizations of the stress field are used for their estimation. 


5.4 Influence of the macroscopic strain level 

In this section the evolution of the parameters of the fitted periodograms as a function of the macro¬ 
scopic strain is investigated. For this purpose the methodology presented in Section [5?^ is applied using 
the realizations of the maximal principal stress fields corresponding to various levels of the loading 
curve, i.e. various values of the equivalent macroscopic strain Eyy = 0., ... , 3.5%. 

The evolution of the standard deviations (cri,cr 2 ) is shown in Figure 17 The two components of 
the periodogram {e.g. Gaussian and exponential) contribute for approximately the same proportion 
to the total variance of the field since the curves are almost superimposed. Note that these standard 
deviations increase with the applied load in the same way as the mean load curve (Figure]^. 

The evolution of the correlation lengths l(x,y)(i, 2 ) is shown in Figure]^ The evolution of the initial 
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Evolution of fx_0 w.r.t the data number 



Evolution of fy_0 w.r.t the data number 



Figure 16: Case #1: Evolution of the fitted initial frequency in the X,Y directions with respect to 
the number of realizations K = 8, ... ,35 


Evolution of sigma w.r.t the macroscopic strain level 



Figure 17: Case #1: Evolution of the fitted standard deviations with respect to the load level (macro¬ 
scopic strain Eyy = 0., ... ,3.5%) 


(1 2 ) 

frequencies is shown in Figure 19 It is observed that once plasticity is settled {i.e. once the 

macroscopic strain Eyy is greater than ^ 0.5%) the parameters describing the fluctuations of the 
maximal principal stress field are almost constant. This conclusion is valid for both the correlation 
lengths and the initial frequencies. Note that the convergence is faster for the parameters related to 
the X direction, i.e. the direction that is transverse to the one-dimensional loading. Finally it is also 
observed that fy^^ is almost equal to zero whatever the load level, thus the zero value in Table 
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Evolution of lx w.r.t the macroscopic strain level 


Evolution of ly w.r.t the macroscopic strain level 




Figure 18: Case #1: Evolution of the fitted correlation lengths in the X, Y directions with respect to 
the load level (macroscopic strain Eyy = 0., ... ,3.5%) 


Evolution of fx_0 w.r.t the macroscopic strain level Evolution of fy_0 w.r.t the macroscopic strain level 




Figure 19: Case #1: Evolution of the fitted initial frequency in the X,Y directions with respect to 
the load level (macroscopic strain Eyy = 0., ... ,3.5%) 
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6 Results — Case #2: random grain geometry 

In this section both the randomness in the grain geometry and in the crystallographic orientations 
are taken into account. A total number of 35 finite element models are run. In each case, the grain 
geometry is obtained from a uniform sampling of points from which a Voronoi tessellation is built. 

6.1 Check of the homogeneity 

As in Section the homogeneity of the maximal principal stress field is checked using the methodology 
proposed in Section [4^ Figure shows the evolution of and CV^j^. These quantities regularly 

decrease and it is seen that they would tend to zero if a larger number of realizations was available. 
This leads to accepting the assumption that the random field is homogeneous. 


Evolution of CV of the empirical mean Evolution of CV of the empirical SD field 




Figure 20: Case #2: Evolution of and CVj^ with respect to the number of realizations 


6.2 Identification of periodograms at 3.5% macroscopic strain 

The average empirical periodogram obtained from L = 35 realizations of the maximal principal stress 
field CT/ at 3.5% of macroscopic strain is plotted in Figure[2^a. Three types of theoretical periodograms 
have been fitted as in the previous section, which lead to the conclusion that the mixed model that 
combines a Gaussian and an exponential component is best suited. The fitted parameters are gathered 
in Table where the parameters fitted for Case #1 are also recalled for the sake of comparison. 

In order to check the accuracy of the fitting, two-dimensional cuts of the empirical (resp. fitted 
periodogram) are plotted in Figure[^(cut along the X direction), Figure[^(cut along the Y direction), 
Figure]^ (cut along the diagonal fx = fy )■ Again the fitting is remarkably accurate, meaning that 
the fitted model of periodogram accurately represents the spatial variability of the maximal principal 
stress field. 

It can be observed from the figures in Table that the fitting is of equal quality in both cases 
(relative error less than 210“^). As far as the contribution of each component of the periodogram is 
concerned, the symmetry reported in Case is not existing anymore since the standard deviation of 
the exponential contribution ((T 2 = 81.6) is much greater than that of the Gaussian part (cti = 35.8). 
The total variance of the field is 7940 MPA^, corresponding to a standard deviation of 89.1 MPa and a 
coefficient of variation of 12%. Thus there is a little more scattering in the random stress field obtained 
in Case ^2 when considering both the random grain geometry and orientations. 

The correlation lengths associated with the exponential part do not differ much in Case #2 com¬ 
pared to Case #1 (corresponding here to 1.5 to 2.4 Dg). In contrast the correlation lengths related 


20 





Figure 21: Case #2: (a) Average empirical periodogram of the stress field at 3.5% macroscopic strain 
- (b) best fitted periodogram 




Figure 22: Case #2: Cut of the periodograms in the X direction 


Table 5: Fitted parameters and error estimates for the mixed “Gaussian + exponential” periodogram 


Case 

e 





Gaussian 




Exponential 


(Eq.( 

32 

0 

CTl 

^xl 

lyl 

AT) 

JxO 

f (J-J 
JyO 

0-2 

1x2 

to 

o to 

JyO 

Case #1 

0.0017 

54.7 

138.4 

159.1 

0.00244 

0 

57.6 

57.5 

63.5 0.00562 

0.0028 

Case #2 

0.0018 

35.8 

269.5 

174.5 

0.00172 

0 

81.6 

67.2 

70.4 0.004 

0 


to the Gaussian part are increased, which tends to produce less rapidly varying realizations. This 
may be explained by the fact that the grain boundaries are “averaged” in Case #2 whereas they were 
fixed in Case #1. The stress concentrations that are usually observed at the grain boundaries are thus 
smoothed in Case #2 compared to Case #1. 
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Figure 23: Case #2: Cut of the periodograms in the Y direction 



• empirical periodogram 
— fitted periodogram 


0.04 0.06 0.08 0.10 


Figure 24: Case #2: Cut of the periodograms along the diagonal = fy 


6.3 Influence of the number of realizations 


In this section one considers the stability of the fitted parameters as a function of the number of 
available realizations K used in the average periodogram method. The procedure applied in the 
previous paragraph is run using if = 8,9, ... , 35 realizations of the stress fiel d. T he evolution of the 
standard deviations (cti, CT 2 ) and the initial frequencies /} 1 ' is shown in Figure 25 (note that /) ( = 0 


in the present case). The evolution of the correlation lengths l(x,y)(i, 2 ) is shown in Figure 26 


From these figures it clearly appears that the fitted parameters are almost constant when the 
number of realizations of the stress field used in their estimation increases. The minimal number of 
K = 8 could be used here without significant errors although it is recommended to keep a value of 
if = 20 as in Case #1 for robustness. 


6.4 Influence of the macroscopic strain level 

Finally the evolution of the parameters of the fitted periodograms as a function of the macroscopic 
strain Eyy is investigated. For this purpose the identihcation method is applied using the realizations 
of the maximal principal stress fields corresponding to various levels of the loading curve, i.e. various 
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Evolution of sigma w.r.t the data number 



Evolution of fx 0 w.r.t the data number 



( 12 ) 

Figure 25: Case #2: Evolution of the fitted standard deviations and the initial frequencies with 
respect to the number of realizations if = 8, ... ,35 

Evolution of lx w.r.t the data number Evolution of ly w.r.t the data number 




Figure 26: Case #2: Evolution of the fitted correlation lengths in the X (left) (resp. Y (right)) 
directions with respect to the number of realizations K = 8, ... ,35 


values of the equivalent macroscopic strain Eyy = 0., ... , 3.5%. 

The evolution of the two standard deviations look similar to the results obtained in Case 
(Figure 27|. It is observed that the ratio uilo\ is almost constant all along the loading path up to 
3.5% strain. As far as the initial frequencies are concerned, there is a complete independance with the 
load level as soon as Eyy is greater than ^ 0.5%, i.e. when plasticity has settled in the aggregate. 
The same conclusion can be drawn for the various correlation lengths. 
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Evolution of sigma w.r.t the macroscopic strain level Evolution of fx_0 w.r.t the macroscopic strain level 



Eyy Eyy 


Figure 27: Case #2: Evolution of the fitted standard deviations (left) with respect to the load level 
(macroscopic strain Eyy = 0., ... ,3.5%) (resp. the initial frequencies (right)) 


Evolution of lx w.r.t the macroscopic strain level Evolution of ly w.r.t the macroscopic strain level 




Figure 28: Case #2: Evolution of the fitted correlation lengths in the X (left) (resp. Y (right)) 
directions with respect to the load level (macroscopic strain Eyy = 0., ... , 3.5%) 
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7 Conclusions 


The distribution of stresses in a material at a microscopic scale (where heterogeneities such as grain 
structures are taken into account) has been given much attention in the context of computational 
homogenization methods. However the current methods usually stick to a deterministic formulation. 
Starting from the premise that any representative volume element (such as a polycrystalline aggregate) 
is a single specific realization of a random quantity, the present paper aims at using methods of 
computational stochastic mechanics for representing the (random) stress field. 

After recalling the basic mathematics of Gaussian random helds, the paper presents a periodogram 
method for estimating the parameters describing the spatial fluctuation of a random field from a 
collection of realizations of this held. This method is adapted in two dimensions from well-known 
techniques originating from signal processing. 

The material under consideration, namely the 16MND5 steel used in nuclear pressure vessels is 
then presented together with a local modelling by polycrystalline hnite element calculations. From 
a collection of 35 realizations of the (maximal principal) stress held, the spatial correlation structure 
of the latter is identihed. By htting various theoretical periodograms, a mixed model combining a 
Gaussian and an exponential-type contribution is retained. These two contributions may be empirically 
interpreted as follows: The Gaussian part corresponds to the huctuation from grain to grain ; the (less 
smooth) exponential component corresponds to the sharp grain boundaries stress concentrations. 

Two cases are considered, namely a “hxed-geometry” case in which only the crystallographic ori¬ 
entations changes within the 35 realizations (hxed grain boundaries), and a “variable geometry” in 
which the grain geometry is randomly sampled for each realization. In both cases, a good convergence 
of the procedure is observed when the number of realizations increases. A set of 20 realizations is 
recommended, although good results are already obtained for ~8 realizations in Gase #2. 

Moreover it is shown that the correlation lengths are of the same order of magnitude as the grain 
size. The initial frequencies that are required for a best htting of the periodogram and that translate 
into some kind of spatial periodicity in the covariogram could be explained by spurious edge effects 
due to the limited size of the aggregate. This should be investigated more in details in further analysis. 

Another important result is drawn from the comparison of the htted parameters at various load 
levels. Once plasticity is settled within the aggregate, the parameters describing the spatial huctuations 
of the held are almost constant. Moreover the variance of the held (sum of the variance of each 
component of the periodogram) increases proportionally to the mean strain/stress curve, meaning that 
the coefficient of variation of the stress held is almost constant (around 11% for the hxed geometry 
and 12% for the variable geometry). 

The results presented in this paper should be conhrmed by additional investigations under different 
types of loading {e.g. biaxial loading). The tools that are presented here may be applicable to three- 
dimensional aggregates and stress helds at a much larger computational cost though. This work is 
currently in progress. 

The identihed stress helds may eventually be re-simulated: new realizations of the stress helds are 
straightforwardly obtained at a low computational cost by random held simulation techniques such as 
the spectral approach or the circulant embedding method (Preumont 1990 Dang et al |2011a|b|. Thi s 
allow us to apply local approach to fracture analysis (such as that presented in Mathieu et al ( 2006| )) 
for the assessment of the brittle fracture of metallic materials, as shown in Dang et al (2013). 
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